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We present a perturbative reconstruction method to make a skymap of gravitational-wave back- 
grounds (GWBs) observed via space-based interferometer. In the presence of anisotropies in GWBs, 
the cross-correlated signals of observed GWBs are inherently time-dependent due to the non- 
stationarity of the gravitational-wave detector. Since the cross-correlated signal is obtained through 
an all-sky integral of primary signals convolving with the antenna pattern function of gravitational- 
wave detectors, the non-stationarity of cross-correlated signals, together with full knowledge of 
antenna pattern functions, can be used to reconstruct an intensity map of the GWBs. Here, we 
give two simple methods to reconstruct a skymap of GWBs based on the perturbative expansion 
in low-frequency regime. The first one is based on harmonic-Fourier representation of data streams 
and the second is based on "direct" time-series data. The latter method enables us to create a 
skymap in a direct manner. The reconstruction technique is demonstrated for the case of the 
Galactic gravitational wave background observed via planned space interferometer, LISA. Although 
' the angular resolution of low-frequency skymap is rather restricted, the methodology presented here 

would be helpful in discriminating the GWBs of galactic origins by those of the extragalactic and/or 
cosmological origins. 



PACS numbers: 04.30.-w, 04.80.Nn, 95.55.Ym, 95.30.Sf 



I. INTRODUCTION 



o 

The gravitational- wave background, incoherent superposition of gravitational waves coming from many unresolved 
\ point-sources and/or diffuse-sources would be a cosmological gold mine to probe the dark side of the Universe. While 
O"^ these signals are generally random and act as confusion noises with respect to a periodic gravitational-wave signal, the 
statistical properties of gravitational- wave backgrounds (GWBs) carry valuable information such as physical processes 
6J[)i of gravitational-wave emission, source distribution and populations. Moreover, the extremely early universe beyond 
the last-scattering surface of cosmic microwave background (CMB) can be directly explored using GWB. Therefore, 
GWBs may be regarded as an ultimate cosmological tool alternative to the CMB. 

Currently, several grand-based gravitational-wave detectors are now under scientific operation, and the search for 
gravitational waves enters a new era. There are several kinds of target GW sources in these detectors, such as, 
coalescence of neutron star binaries and core-collapsed supernova. On the other hand, planned space-based detector, 
LISA (Laser Interferometer Space Antenna) and the next-generation detectors, e.g., DECIGO [l] and BBO aim 
at detecting gravitational waves in low-frequency band O.lmHz -0.1Hz, in which many detectable candidates for 
GWBs exist. Among them, the GWB originated from the inflationary epoch may be detected directly in this band 
(e.g., 0,0|). Therefore, for future application to cosmology, various implications for these GWBs should deserve 
consideration both from the theoretical and the observational viewpoint. 

One fundamental issue to access a new subject of cosmology is to make a skymap of GWB as a first step. Similar 
to the case of the CMB, the intensity map of the GWBs is of particular interest and it provides valuable information, 
which plays a key role to clarify the origin of GWBs. Importantly, the low-frequency GWBs observed by LISA are 
expected to be anisotropic due to the contribution of Galactic binaries to the confusion noise (0, 0, 0i see a l so 
recent numerical simulations H,|UEi)- Thus, the GWB skymap is potentially useful to discriminate the individual 
backgrounds from many superposed GWBs, as well as to identify the underlying physical processes. The basic idea to 
create a skymap of GWB is to use the directional information obtained through the time-modulation of the correlation 
signals, which is caused by the motion of gravitational-wave detector. As detector's antenna pattern sweeps out the 
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sky, the amplitude of the gravitational-wave signal would gradually change in time if the intensity distribution of 
GWB is anisotropic. Using this, a method to explore anisotropies of GWB has been proposed [TT1 IT3 . IT3 . Later, the 
methodology was applied to study the anisotropic GWB observed by space interferometer, LISA [l4 IWL IWL Il7l llSl] . 

In the previous paper which is referred to as Paper I in the present paper, we have investigated the directional 
sensitivity of space interferometer to the anisotropy of GWB. Particularly focusing on the geometric properties of 
antenna pattern functions and their angular power, we found that the angular sensitivity to the anisotropic GWBs 
is severely restricted by the data combination and the symmetry of detector configuration. As a result, in the case of 
the single LISA detector, detectable multipole moments are limited to £ < 8-10 with the effective strain sensitivity 
h ~ 10 -20 Hz -1 / 2 . This is marked contrast to the angular sensitivity to the chirp signals emitted from point sources, 
in which the angular resolution can reach at a level of a square degree or even better than that [23, El H3, • 

Despite the poor resolution of LISA detector with respect to GWBs, making a skymap of GWBs is still an important 
general issue and thus needs to be investigated. In the present paper, we continue to investigate the map-making 
problem and consider how the intensity map of the GWBs is reconstructed from the time-modulation signals observed 
via space interferometer. In particular, we are interested in the low-frequency GWB, the wavelength of which is 
longer than the arm-length of the detector. In such case, the frequency dependence of the detector response becomes 
simpler and a perturbative scheme based on the low-frequency expansion of the antenna pattern functions can be 
applied. Owing to the least-squares approximation, we present a robust reconstruction method. The methodology is 
quite general and is also applicable to the map-making problem in the case of the ground detectors. We demonstrate 
how the present reconstruction method works well in a specific example of GWB source, i.e., Galactic confusion- noise 
background. With a sufficient high signal-to-noise ratio for each anisotropic components of GWB, we show that the 
space interferometer LISA can create the low- frequency skymap of Galactic GWB with angular resolution I < 5. 
Since the resultant skymap is obtained in a non-parametric way without any assumption of source distributions, we 
hope that despite the poor angular resolution, the present methodology will be helpful to give a tight constraint on 
the luminosity distribution of GWBs if we combine it with the other observational techniques. 

The organization of the paper is as follows. In Sec. |H] a brief discussion on the detection and the signal processing 
of anisotropic GWBs is presented, together with the detector response of space interferometer. Scc lIIII describes the 
details of the reconstruction of a GWB skymap in low-frequency regime. Owing to the least-squares approximation 
and the low-frequency expansion, a perturbative reconstruction scheme is developed. Related to this (in the case 
of LISA), and we give an important remark on the degeneracy between some multipole moments (Appendix [BJ. In 
Sec lIVI the reconstruction method is demonstrated in the case of the Galactic GWB. The signal-to-noise ratios for 
anisotropy of GWB are evaluated and the feasibility to make a skymap of GWB is discussed. Finally, Section Ivl is 
devoted to a summary and conclusion. 



II. BASIC FORMALISM 
A. Correlation analysis 

Let us begin with briefly reviewing the signal processing of gravitational-wave backgrounds based on the corre- 
lation analysis [19|. Stochastic gravitational- wave backgrounds are described by incoherent superposition of plane 
gravitational waves h = hij given by 

/+oo p 
df / dn 2 -/(*-«-) h A (f, n) e A (n). (i) 

A=-t-,X "°° J 

The Fourier amplitude h,A(f, fi) of the gravitational waves for the two polarization modes e A (A = +, x) is assumed 
to be characterized by a stationary random process with zero mean (h A ) = 0. The power spectral density Sh is then 
defined by 

(h* A (f,si)h A ,(f,n')) = x -8{f - f f {n ~ n>) 5AA'S h (\f\, a), (2) 

where $7 is the direction of a propagating plane wave. Note that the statistical independence between two different 
directions in the sky is implicitly assumed in equation J2Jl, which might not be generally correct. Actually, CMB 
skymap exhibits a large-angle correlation between the different skies, which reflects the primordial density fluctuations. 
For the GWB of our interest, the wavelength of the tensor fluctuations detected via space interferometers is much 
shorter than the cosmological scales and thereby the assumption put in equation is practically valid. 

The detection of a gravitational- wave background is achieved through the correlation analysis of two data streams. 
The planned space interferometer, LISA and also the next generation detectors DECIGO/BBO constitute several 
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spacecrafts, each of which exchanges laser beams with the others. Combining these laser pulses, it is possible to 
synthesize the various output streams which are sensitive (or insensitive) to the gravitational-wave signal. The output 
stream for a specific combination / denoted by sj (t) is generally described by a sum of the gravitational- wave signal 
hi(t) and the instrumental noise ni(t) by 

Sl (t) = hi(t) + m(t). 

We assume that the noise ni(t) is treated as a Gaussian random process with spectral density S n (f) and zero mean 
(nj) = 0. The gravitational- wave signal hj is obtained by contracting h with detector's response function and/or 
detector tensor D := D^: 

hi® = / +OC d f I dn e i27r/(t -"- Xj3 D/(n, /; t) : e A (n) h A (f, ft). (3) 

Note that the response function explicitly depends on time. The time variation of response function is caused by the 
orbital motion of the detector and it plays a key role in reconstructing a skymap of the GWBs. Typically, the orbital 
frequency of the detector motion is much lower than the observed frequency and within the case, the expression © 
is validated. 

Provided the two output data sets, the correlation analysis is examined depending on the strategy of data analysis, 
i.e., self-correlation analysis using the single data stream or cross-correlation analysis using the two independent 
data stream. Defining Su(t) = (sj(t)sj(t)), the Fourier counterpart of it Sjj(t,f), which is related with Sjj(t) by 
Su(t) = / dfSu(t,f), becomes 1 

Su(t J) = Cuitjj + Su S n (\f\), (4) 

where we define 

CuitJ) = J ^s h (\f\,n) Tfj{f,n- 1). (5) 

Here, the function J-fj is the so-called antenna pattern function defined in an ecliptic coordinate, which is expressed 
in terms of detector's response function: 



A= + ,x 

F I (nj;t)=r> I (n,f;t):e A (tl) (6) 



Apart from the second term, equation (JIJ implies that the luminosity distribution of GWBs Sh(f, O) can be obtained 
by deconvolving the all-sky integral of antenna pattern function from the time-series data Su(t). To see this more 
explicitly, we focus on equation and decompose the antenna pattern function and the luminosity distribution into 
spherical harmonics in an ecliptic coordinate, i.e., sky-fixed frame: 

s h (\f\,n) = J2\ptm(f)]*Y e * m (si), 
rfjtf, n-,t) = J2 t) Y tm (a). (7) 

Note that the properties of spherical harmonics yield p\ m = (— l)' m pi.- m and a\ m = (— f ) f ~ m a^_ m , where the latter 
property comes from Ffjif, Sl;t) = Fuif, — Q;t) Jjj. Substituting into JSJ becomes 

/) = J 1 £ \P*™(f)T a f™(f, *), (8) 
im 

where we have dropped the contribution from the detector noise. The expression (JS| still involves the time-dependence 
of antenna pattern functions. To eliminate this, one may further rewrite equation © by employing the harmonic 



1 In Paper I, we had used Cij to denote the output data (si(t)sj(t)) itself. This coincides with the present definition if we neglect the 
instrumental noises. 
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FIG. 1: LISA configuration in the ecliptic coordinates. The Galactic center is at RA 267.4° and dec. —5.5° in this coordinate 
system. 



expansion in detector's rest frame. We denote the multipole coefficients of the antenna pattern in detector's rest frame 
by ai m . The transformation between the detector rest frame and the sky- fixed frame is described by a rotation matrix 
by the Euler angles (V>, <p), whose explicit relation is expressed in terms of the Wigner D matrices 2 [TliL Il6l l24| : 

t 

af m (f,t) = £ e-^di^e-^ a tn (f). (9) 

Here the Euler rotation is defined to perform a sequence of rotation, starting with a rotation by ip about the original 
z axis, followed by rotation by •& about the original y axis, and ending with a rotation by if about the original z axis. 
The Wigner D matrices for n > to is 



^(0) = (-!)' "y (£ + w)! ( £ _ m) , p in- H-cose) (10) 

with Pn a ' b ^ being Jacobi polynomial. For n < to, we have d l nm = (— l) n ~ m d^ nn . 

Let us now focus on the orbital motion of the LISA constellation (see Fig^l. The LISA orbital motion can be 
expressed by ip = —ut, = — 7r/3, tp — cut, where lo = 2t:/Tq is the orbital frequency of LISA (To = 1 sidereal year) 
3 Since the antenna pattern function is periodic in time due to the orbital motion, the expected signals also vary in 
time with the same period T . It is therefore convenient to express the output signals by 

Cu(t,f)= CuAf)e* kut . (11) 

k— — oo 

Using the relation with specific parameter set, the Fourier component Cjj^if) is then given by 

Ck(f) = i [ T ° 'dte~ ikut C u (t,f) 

^ x i— k 

= i^Yl Yl \Pe™(f)}* d\m+k),m (~|) a L{m+k)U) (12) 



=0 m=-t 



for k > 0. As for k < 0, the lower and the upper limit of the sum over to are changed to to = — i — k and to = £, 
respectively. 



2 We are specifically concerned with the time dependence of directional sensitivity of antenna pattern functions, not the real orbital 
motion. Hence, only the time evolution of directional dependence is considered in the expression 

3 The relation tp = — ip does not necessarily hold for orbital motion of LISA and there may be some possibilities to impose a constant 
phase difference, i.e., ip = —if + c. However, for the sake of simplicity, we put c = 0. 
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Equation l(T^)l as well as (J5J is the theoretical basis to reconstruct the skymap of GWBs. Given the output data 
Ck{f) (or Cij(t, /)) experimentally, the task is to solve the linear system l|12[) with respect to pe m (f) for given antenna 
pattern functions. One important remark deduced from equation l|12|l is that the accessible multipole coefficients pi m 
are severely restricted by the angular sensitivity of antenna pattern functions. The important properties of the antenna 
pattern functions are summarized in next subsection. Another important message is that the above linear systems 
are generally either over-constrained or under-determined. For the expression 1)12(1 . if one truncates the multipole 
expansion with I = £ max , the system consists of \{t m ax + l){^max + 1) unknowns and N(2£ max + 1) equations (for 
k > 0), where N is the number of available modes of the antenna pattern functions. Thus this deconvolution problem 
is, in principle, over-determined for a relatively small truncation multipole ^ max , while it becomes under-determined 
for a larger value of £ max . We will later discuss the over-determined case in Sec II V C"T1 where £ max = 5 and N ~ 5. 



B. Detector response and antenna pattern functions 



The output signals of space interferometer sensitive to the gravitational waves are constructed by time-delayed 
combination of laser pulses. In the case of LISA, the technique to synthesize data streams canceling the laser 
frequency noise is known as time-delay interferometry (TDI), which is crucial for our subsequent analysis. In the 
present paper, we use the optimal set of TDI variables independently found by Prince et al. [23 and Nayak et al.[26|. 
which are free from the noise correlation 4 . A simple realization of such data set is obtained from a combination of 
Sagnac signals, which are the six-link observables using all six LISA oriented arms. For example, the Sagnac signal Si 
measures the phase difference accumulated by two laser beams received at space craft 1, each of which travels around 
the LISA array in clockwise or counter-clockwise direction. The explicit expression of detector tensor for such signal, 
which we denote by D Sl , is given in equation (21) of Paper I (see also [23, The analytic expression for other 

detector tensors Dg 2 and Ds :i are also obtained by the cyclic permutation of the unit vectors a, b and c. 

Combining the three Sagnac signals, a set of optimal data combinations can be constructed [25|,|26| (see also (2!j): 

D A =-^(D S3 -D Sl ), 

D E = i=(Ds 1 -2Ds 2 +Ds3), 

D T = -^=(D Sl +D S2 +D S3 ). (13) 

These three detector tensors are referred to as A, E, T-variables, which generate six kinds of antenna pattern functions. 
Notice that in the equal arm- length limit, frequency dependence of these functions is simply expressed in terms of the 
normalized frequency: 

(») 

where the characteristic frequency is given by /» = c/(2irL) with L being the arm-length of detector. With the 
arm-length L = 5x 10 6 km, the characteristic frequency of LISA becomes /» ~ 9.54 mHz. In what follows, we use the 
analytic expressions for equal arm-length limit to demonstrate the reconstruction of GWB skymap. This is sufficient 
for the present purpose, because we are concerned with a fundamental theoretical basis to map-making capability of 



TABLE I: Important properties of antenna pattern function. Note that the cross-correlated data are blind to isotropic GWBs. 



Combination of variables Visible multipole moments in low-frequency regime (/ < 1) General properties 



(A, A), (E,E) 


O(P) : 
0(f 4 ) : 
0(f 2 ) : 


1 = 0,2, 


4 


visible only to I =even 


(T,T) 


e = o, 4, 


6 


visible only to I =even 


(A,E) 


1 = 4, 


0(/ 3 ) : £ = 3, 5 


blind to I = 0, 1 


(A,T), (E,T) 


0(f 3 ) : 


1 = 1,3, 


5 


blind to I = 



4 The optimal TDI variables adopted here may not be the best TDI combinations for the present purpose. There might be a better choice 
of the signal combinations, although a dramatic improvement of the angular sensitivity would not be expected. 
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GWBs. The idea provided in the present paper would allow us to examine a more realistic situation and the same 
strategy can be applied to an extended analysis in the same manner. 

Finally, we note that the antenna pattern functions for the optimal combinations of TDI have several distinctive 
features in angular sensitivity. In the low frequency limit (/ <C 1), the antenna pattern functions for the self- 
correlations and the cross-correlations are expanded as 



T AA (f, n)= ^imf 2 




+ ^Ami A 


+ o(/ 6 


TttU, n) = 




4t(^)/ 4 


+ o(/ 6 


T AE (f,n) = T A 2) E (si)p 


+■ ^llmP 


+ ^Ito/ 4 


+ o(f 


T AT {.f, O) = 






+ 0(f 5 



(15) 



The frequency dependence of Tee and Tet are the same as for Taa and Tat, respectively. Since the leading term of 
the TT-correlation is 0(/ 4 ), the TT-correlation becomes insensitive to the gravitational waves in the low-frequency 
regime. We will use all correlated data set except for the TT-correlation. In Table we summarize the important 
properties for the multipole moments of antenna pattern functions. Also, in Appendix lAl employing the perturbative 
approach based on the low- frequency approximation /«1, the spherical harmonic expansion for antenna pattern 
functions are analytically computed, which will be useful in subsequent analysis. 

III. PERTURBATIVE RECONSTRUCTION METHOD FOR GWB SKYMAP 

We are in position to discuss the methodology to reconstruct a skymap of GWB based on the expression l|12|) 
(or (JSJ). Since we are specifically concerned with low- frequency sources observed via LISA, it would be helpful to 
employ a perturbative approach using the low-frequency expansion of antenna pattern function. In Sec lIII Al owing to 
the expression (|12f) . a perturbative reconstruction method is presented. Sec lIII Bl discusses alternative reconstruction 
method based on the time-series representation ||5J. 

A. General scheme 

Hereafter, we focus on the GWBs in the low-frequency band of the detector, the wavelength of which is typically 
longer than the arm- length of the gravitational detector, i.e., / < 1. Without loss of generality, we restrict our 
attention to the reconstruction of a GWB skymap in a certain narrow frequency range, / ~ / + A/, within which a 
separable form of the GWB spectrum is a good assumption: 

S h (f,Sl)=H(f)P(n). (16) 

Further, for our interest of the narrow bandwidth, it is reasonable to assume that the spectral density H(f) is described 
by a power-law form as H(f) = M f a . In the reconstruction analysis discussed below, the spectral index a is assumed 
to be determined beforehand from theoretical prediction and/or experimental constraint 5 . 

In the low frequency approximation up to the order C(/ 3 ), we have five output signals which respond to the GWBs, 
i.e., AA-, EE-, AE-, AT- and £T-corrclations, each output of which is represented by equation l|12|). Collecting these 
linear equations and arraying them appropriately, the linear algebraic equations can be symbolically written in a 
matrix form as 

c(/)=A(/)-p. (17) 

In the above expression, while the vector p represents the unknowns consisting of the multipole coefficients of GWBs 
•pim, the vector c(/) contains the correlation signals Ck(f)- Here, the frequency dependence of the function H(f) has 
been already factorized and thereby the vector p only contains the information about the angular distribution. The 



5 In our general scheme, we do not assume a priori the normalization factor J\[ in the function H(f), which should be simultaneously 
determined with the reconstruction of an intensity distribution P(fl). In the low- frequency approximation, / < 1, however, there exists 
a degeneracy between the monopole and the quadrupole components and one cannot correctly determine the normalization M ■ (See 
Appendix iBll 
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matrix A(/) is the known quantity consisting of the multipole coefficients of each antenna pattern and the Wigner 
D matrices (see Appendix IO for an explicit example). 

As we have already mentioned, the linear systems (|17f) become either over-determined or under-determined system. 
In most of our treatment in the low-frequency regime, the linear systems (|17|l tend to become over-determined, but 
this is not always correct depending on the amplitude of GWB spectrum H(f) (see Sec lIV Bj) . In any case, the matrix 
A would not be a square matrix and the number of components of the vector c does not coincide with the one in the 
vector p. In the over-determined case, while there is a hope to get a unique solution p to produce the gravitational- 
wave signal c, it seems practically difficult due to the errors associated with the instrumental noise and/or numerical 
analysis. Hence, instead of pursuit of a rigorous solution, it would be better to focus on the issue how to get an 
approximate solution by a simple and systematic method. 

In the case of our linear system, the approximate solution for the multipole coefficient p apP rox can be obtained from 
the least-squares method in the following form: 

p_ - A+ ■ c. (18) 

The matrix A + is called the pseudo-inverse matrix of Moore-Penrose type, whose explicit expression can be uniquely 
determined from the singular value decomposition (SVD) [^J. According to the theorem of linear algebra, the 
matrix A, whose number of rows is greater than or equal to its number of columns, can be generally written as 
A = U T ■ diag[wi] • V . Here, the matrices U and V are orthonormal matrices which satisfy W U = V' V = 1, where 
the quantity with subscript t represents the Hermite conjugate variable. The quantity diag[u>i] represents an diagonal 
matrix with singular values tUj with respect to the matrix A. Then the pseudo-inverse matrix becomes 

A+ = V T ■ diagK -1 ] • U. (19) 

It should be stressed that the explicit form of the pseudo-inverse matrix A + is characterized only by the angular 
dependence of antenna pattern functions. 

In principle, the least-squares method by SVD can work well and all the accessible multipole moments of GWB 
would be obtained as long as the antenna pattern functions have the corresponding sensitivity to each detectable 
multipole moment. As we mentioned in Sec III Bl however, the angular power of antenna pattern function depends 
sensitively on the frequency. In the low-frequency regime, the frequency dependence of the non-vanishing multipole 
moments appears at 0(f 2 ) for I = 0, 2 and 4, and 0(f 3 ) for I = 1, 3 and 5 (see Table QJ. In this respect, by a 
naive application of the least-squares method, it is difficult to extract the information about £ =odd modes of GWBs 
because the singular values of the matrix A arc dominated by the lowest-order contribution of the antenna pattern 
functions. 

For a practical and a reliable estimate of the odd multipoles in low-frequency regime, the least-squares method 
should be applied combining with the perturbative scheme described below. Let us recall from the expression i|15|) 
that the matrix A can be expanded as 

A = f 2 + /* A^ + • • • (20) 
for a matrix consisting of the self-correlation signals T aa and Tee, 

A = f 2 A^ + f 3 AW + / 4 A (4) + . . . ( 2i) 
for a matrix consisting of the cross-correlation signal Tae, and 

A = f 3 A^ + ,/ 4 A< 4 > + • • • (22) 

for a matrix consisting of the cross-correlation signals Tat and Tet- Then the resultant matrices become 
independent of the frequency. The above perturbative expansion implies that the output signal c is also expanded in 
powers of /. Since the frequency dependence of the function H(f) has been already factorized in equation (|17|) . we 
have 

f 2 + / 4 cW + • • • , for AA-, EE-correlations 

c(/) = I f 2 c& + f 3 c< 3 ) + /V 4 ) + • • • , for AE-correlation (23) 
f 3 c (3 > + /V 4 ) + • • • , for AT-,ET-correlations 
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Substituting the terms 120l) - (|22|l and ill'MI) into the expression l|17|l and collecting the terms of each order of /, we have 



(i = 2, 3, 4, 



(24) 



where the subscript W means the quantity consisting of the order 0(f l ) terms. The vector pW represents the accessible 
multipolc moments pg m to which the antenna pattern functions become sensitive in this order. For example, the vector 
p' 2 ' contains the I = 0, 2 and 4 modes, while the vector p^ 3 ^ have the multipole moments with £ = 1, 3 and 5. Since 
the expression (|24|l has no explicit frequency dependence, we can reliably apply the least-squares solution by the SVD 
to reconstruct the odd modes of GWBs, as well as the even modes: 



(25) 



Each step of the perturbative reconstruction scheme is summarized in Fig. [2] (see Appendix IO in more details). 



multipole moments of antenna pattern 
functions at detector's rest frame 



aim(f) 



expand in powers of f 



ain 



A 



(i) 



pseudo-inverse matrix 



(see Appendix A) 



5 output data 
(AA, EE, AE, AT, ET) 



Wigner D matrices 

df m (0) 



Fourier components 
of correlation signals 



/multipole coefficientsN. Tp^l ppIQX zz ] ' C 

of GWB skymaj^ 

multiply the function H(f) 
finally 



C k (f) 



expand in powers of f 



Lk 



factorize the frequency 

dependence of H(f) 



FIG. 2: Flowchart of perturbative reconstruction scheme based on the harmonic-Fourier representation. 

Finally, we note that the low-frequency reconstruction method presented here assumes the perturbative expansion 
form of the output signal c. To determine the coefficients cW in equation (|23|l . one must know the frequency 
dependence of the vector c in the narrow bandwidth / ~ / + A/, which can be achieved by analyzing the multi- 
frequency data. One important remark is that the signal-to-noise ratio in each reconstructed multipoles might be 
influenced by the sampling frequencies and/or the data analysis strategy. This point will be discussed later in Sec. 

ITvcl 



B. Reconstruction based on the time-series representation 



So far, we have discussed the reconstruction method based on the harmonic-Fourier representation (|12f> . The main 
advantage of the harmonic-Fourier representation is that it gives a simple algebraic equation suitable for applying 
the least-squares solution by SVD. However, the harmonic-Fourier representation implicitly assumes that the space 
interferometer orbits around the sun under keeping their configuration rigidly. In reality, rigid adiabatic treatment of 
the spacecraft motion is no longer valid due to the intrinsic variation of arm-length caused by the Keplerian motion 
of three space crafts, as well as the tidal perturbation by the gravitational force of solar system planets [3ll l32l l33| . 
In a rigorous sense, the time dependence of antenna pattern function cannot be described by the Euler rotation of 
antenna pattern function at rest frame (see footnote [2J . Although the influence of arm- length variation is expected 
to be small in the low-frequency regime, alternative approach based on the other representation would be helpful to 
consider the GWB skymap beyond the low-frequency regime. 
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Here, we briefly discuss the reconstruction method based on the time-series representation The time-series 
representation is mathematically equivalent to the harmonic-Fourier representation under the rigid adiabatic treat- 
ment, but in general, no additional assumption for space craft configuration is invoked to derive the expression 
except for the premise that the time-dependence of the antenna pattern functions, i.e., the motion of each spacecraft, 
is well-known theoretically and/or observationally. Moreover, as we see below, the method based on the time-series 
representation allows us to reconstruct directly the skymap Sh, without going through intermediate variables, e.g. 
Pi m . In this sense, the time-series representation would be superior to the harmonic- Fourier representation for prac- 
tical purpose, although there still remains the same problem as discussed in Appendix [B] concerning the degeneracy 
of the multipole coefficients. 

In principle, the perturbative reconstruction scheme given in Sec llll Al can be applicable to the map-making problem 
based on the time-series representation. To apply this, we discretize the integral expression © so as to reduce it to 
a matrix form like equation l|17l) . For example, we discretize the celestial sphere (8, </>) into a regular TV meshes and 
also the continuous time-series into the regular M grids. Then we have 



N AO 

c{t u f) = V s h (\f\, n,) T E (f, n i5 u) — >- 



47T 



(t = l,2,...,M), 



(26) 



where AO^ = sin 9i AOi A</>j. Here, we have ignored the noise contribution and dropped the subscript jj. For a 
reconstruction of low-frequency skymap, a large number of mesh and/or grid are not necessary. Typically, for the 
angular resolution up to £ = 5, it is sufficient to set M = 16 and N = 16 (see Sec lIV C|l . The above discretization 
procedure is repeated for the five output data of the correlation signals. Then, simply following the same procedure as 
in the case of the harmonic-Fourier representation, the least-squares method by SVD is applied to get the luminosity 
distribution of GWBs. Note that the meanings of the matrix A^, the vectors c'*' and in the perturbative 
expansion are appropriately changed as follows. For the matrix we have 



A« = — 

47T 



V (tit, t M ) a«i (fh,t M ) An 2 



T^(n N ,t 2 )An N 
^(n N ,t M )An N J 



(4 = 2,3, •••) 



(27) 



The corresponding vectors C(j) and p(j) becomes 

cW(ta) 



»w = 



( s^ono \ 
^°(n 2 ) 



(28) 



where C'^(t) and S^(Cl) the perturbative coefficients of C(t,f) and S/j(|/|,f2) in power of p, respectively. Using 

(2) (31 

these expressions, the least-squares solution is constructed as p apP rox = Papprox + Papprox -I , with a help of equation 

(J2U. Then, the resultant expression p apP rox directly gives a GWB skymap in the ecliptic coordinate, i.e., Sh(\f\, 
not the multipole coefficients. 



IV. DEMONSTRATION: SKYMAP OF GALACTIC BACKGROUND 



Perturbative reconstruction scheme presented in the previous section is applicable to a general map-making problem 
for any kind of GWB sources. In this section, to see how our general scheme works in practice, we demonstrate the 
reconstruction of a GWB skymap focusing on a specific source of anisotropic GWB. An interesting example for LISA 
detector is a confusion-noise background produced by the Galactic population of unresolved binaries. After describing 
a model of Galactic GWB in Sec lIV Al the expected signals for time-modulation data of self- and cross-correlated 
TDIs are calculated in Sec lIV Bl Based on these, the detectable Fourier components for time-modulation signals are 
discussed evaluating the signal-to-noise ratios. In Sec lIV Cl a reconstruction of the GWB skymap is performed based 
on the harmonic-Fourier representation and the time-series representation. The resultant values of the multipole 
coefficients for Galactic GWBs are compared with those from the full-resolution skymap, taking account of the 
influence of the noises. 
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A. A model of Galactic GWB 



For our interest of GWBs in the low-frequency regime with / < /* ~ 9.54 mHz, it is reasonable to assume that the 
anisotropic GWB spectrum Si l (f,tt) is separately treated as Sh(f,Sl) = H(f) P(tt), as we mentioned. The power 
spectral density H(f) is approximated by a power-law function like H(f) = N/ a , whose amplitude will be explicitly 
given later. For illustrative purpose, we consider the simplest model of luminosity distribution -P(J~2), in which the 
Galactic GWB is described by an incoherent superposition of gravitational waves produced by compact binaries whose 
spatial structure just traces the Galactic stellar distribution observed via infrared photometry. We use the fitting 
model of Galactic stellar distribution given in j3^|, which consists of triaxial bulge and disk components (see also [l7|L 
The explicit functional form of the density distribution p(x) written in the Galactic coordinate system becomes 

p(x) = pdisk(x) + Pbuigc(x) ; (29) 



Pbulgc ~ (1 + a/ao) 1 - 8 

e -\z\/z g— \z\/zi 



Pdisk = r a 

z Z\ 



with the quantities a and R defined by a = [x 2 + {y/ri) 2 + (z/C) 2 ] 1 ^ 2 and R = (x 2 + y 2 ) 1 ^ 2 - Note that the z-axis 
is oriented to the north Galactic pole and the direction of the rr-axis is 20° different from the Sun-center line. Here 
the parameters are given as follows: po = 624, a m = 1.9 kpc, ao = 100 pc, R s = 2.5 kpc, zq = 210 pc, z\ — 42 pc, 
a = 0.27, r\ = 0.5 and C = 0.6. Provided the three-dimensional structure of stellar distribution p(x), the angular 
distribution -P(Jl) is obtained by projecting it onto the sphere in observed frame, i.e., ecliptic coordinate: 

P(fl)=C / drAnr 2 (30) 

where C is a numerical constant normalized by J dft P(ft) = 1. The integration over r is performed along a line-of- 
sight direction from a location of space interferometer to infinity. 

In FigOU the projected intensity distribution of GWB is numerically obtained specifically in ecliptic coordinate and 
it is then transformed into Galactic coordinate, shown as the Hammer-Aitoff map. A strong intensity peak is found 
around the Galactic center and the disk-like structure can be clearly seen. Notice that while the result depicted 
in Fig|31 represents a full-resolution skymap, it cannot be attained from the perturbative reconstruction scheme in 
low-frequency regime. For the antenna pattern functions of cross-correlated TDI signals up to the order C(/ 3 ), the 
detectable multipole moments of GWB are limited to £ < 5. Moreover, in the low-frequency limit 0(f 2 ), only the 
even modes £ — 0, 2 and 4 are measurable. Thus, the reconstructed skymap would be rather miserable compared 
to the full skymap which contains the higher multipoles I > 30 (see FigEl m Appendix). Taking account of these 
facts, in right panel of Fig^l we plot the low-resolution skymap, which was obtained from the full-resolution skymap 
just dropping the higher multipole moments with I > 5. In Appendix [Dj with a help of the Fortran package of 
spherical harmonic analysis ([35|], see Appendix ITjl). the numerical values of the multipole coefficients pe m up to I = 5 
are computed and are summarized in Table IIHI Also in the left panel, the odd modes are further subtracted and 
the remaining multipole moments are only £ = 0, 2 and 4. As a result, the fine structure around the bulge and the 
disk components is coarse-grained and the intensity of the GWB diminishes. It also shows some fake patterns with 
negative intensity. Nevertheless, one can clearly see the anisotropic structure of GWB , which is mainly contributed 
from the gravitational-wave sources around the Galactic disk. With a perturbative reconstruction of low-frequency 
up to C(/ 3 ), one can roughly infer that the main sources of Galactic GWB comes from the Galactic center. 



B. Time-modulation signals and signal-to-noise ratio 

Given the intensity distribution of GWB, one can calculate the cross-correlation signals observed via LISA, which 
are inherently time-dependent due to the orbital motion of the LISA detector. The effect of the annual modulation of 
the Galactic binary confusion noise on the LISA data analysis had been previously studied in |17| in the low-frequency 
limit /« 1. Recently, Monte Carlo simulations of Galactic GWB were carried out by several groups and the annual 
modulation of GWB intensity has been confirmed in a realistic setup with specific detector output |3, III [lCj . 

Owing to the expression JSJ, the time-modulation signals C(£, /) neglecting the instrumental noises are computed 
for optimal TDIs at the frequency / = 0.1, i.e., / ~ lmHz using the full expression of antenna pattern function ijfJJl. 
The results are then plotted as function of orbital phase. In Fig. [3 six outputs of the self- and cross-correlation 
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FIG. 3: Full-resolution skymap of the Galactic gravitational-wave background shown as Hammer- Aitoff map in Galactic 
coordinate. 
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FIG. 4: Expected images of GWB skymap by LISA based on the low-frequency reconstruction scheme, which are both depicted 
in the Galactic coordinate. Left: low-resolution skymap created by truncating all the multipoles except for I — 0, 2 and 4 from 
the original full-resolution skymap. Right: skymap created by truncating the higher multipoles £ > 6. 



signals normalized to its time-averaged value, C(t, /)/|C*o(/)| are shown. The solid and dashed lines denote the real 
and imaginary parts of the correlation signals, respectively. 

The time modulation of these outputs basically reflects the spatial structure of GWB. As LISA orbits around 
the Sun, the direction normal to LISA's detector plane, which is the direction sensitive to the gravitational waves, 
sweeps across the Galactic plane. Since the response of the LISA detector to the gravitational waves along the Cl 
direction give the same response to the waves along the —fl direction (see Eq.Q and the brief comment there), the 
time-modulation signal is expected to have a bimodal structure, like AA- and E'E'-correlations. However, the actual 
modulation signals are more complicated, depending on the angular sensitivity of their antenna pattern functions, 
as well as the observed frequency. Further, the cross-correlation data can be generally complex variables, whose 
behaviors are different between real- and imaginary-parts. 

Notice that the time- modulation signals presented in Fig. [3] are the results in an idealistic situation free from the 
noise contribution. In presence of the random noise, some of the correlation data which contain gravitational wave 
signals are overwhelmed by noise, which cannot be used for the reconstruction of GWB skymap. Hence, one must 
consider the signal-to-noise ratio (SNR) to discriminate available correlation data. To evaluate this, the output signals 
depicted in Fig[31are first transformed to its Fourier counterpart, Ck(f)- The resultant Fourier components for each 
signal are then compared to the noise contributions. The SNR for each component is expressed in the following form 
(|l|, see also [HEI): 

(D.-^^f- < 31) 

Here we set the observational time to T = 10 8 sec and the bandwidth to A/ = 10 _3 Hz. The quantity Nk represents 
the noise contribution. The important remark is that the noise contribution comes from not only the detector noise 
but also the randomness of the signal itself. The details of the analytic expression for Nk will be given elsewhere |3G§ . 
Here, as a crude estimate, we evaluate the noise contribution Nk as 



N k = S?, (I = A, E, T) 



(32) 
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FIG. 5: Annual time-modulation signals of self- and cross-correlation data as function of orbital phase assuming the observed 
frequency / = 0.1, i.e., / ~ lmHz. Here, the correlation signals Cij(t,f) are plotted normalizing with k — component of 
correlation signal Cjj,o(/). Solid and dashed lines represent the real and the imaginary part of correlation signal, respectively. 



for k = component of self-correlation signals and 



N k = ^max (c u , Cjj, , C U , Q S£ J , Cjj^Sf/, Spsrf , (I, J = A, E, T) (33) 

for cross-correlation signals and k ^ component of self-correlation signals. To estimate SNR, one further needs the 
spectral density for instrumental noise, i.e., S^ A (f), S EE (f) and S^ T (f). We use the expressions given in equation 
(58) of Paper I (see also USES]) 6 . Note that all the cross-correlated noise spectra such as S£ E (f) and S AT (f) are 
exactly canceled. 

In FigEJ the SNRs for six output data are evaluated and are shown in the histogram as function of Fourier 
component, k. In these panels, thick-dotted lines show the detection limit of (S/N)k = 5, while the thin-dotted lines 
mean (S/N)k = 1. When evaluating the SNR, we specifically consider the two cases: 

Case A: realistic case in which the rms amplitude of GWB spectrum is given by Sf/ 2 = 5 x lfr^Hz- 1 / 2 at the 
frequency / = lmHz 7 . 

Case B: optimistic case in which the rms amplitude of GWB spectrum is ten times larger than that in the realistic 
case, i.e., S 1 , 1 / 2 = 5 x lCT^Hz- 1 / 2 at / = lmHz. 

The results of SNR are then shown in solid (case A) and dashed lines (case B), respectively. 

As anticipated from the sensitivity of antenna pattern function in Tabled the SNR for TT-correlation is much less 
than unity. Even in the optimistic case, the SNR is about ten times smaller than unity and thus the TT-correlation 
data cannot be used for reconstruction of the GWB skymap. Apart from this, the SNRs for the self-correlation signals 



6 In equation (58) of Paper I, there are some typos in the numerical values of S s h Q t(f) and 5 arrf ,i (/ ). In the present paper, we adopt 
S'shotC/) = 16 x 10~ 41 Hz- 1 and 5 acccl (/) = 2.31 X 1CT 41 (mHz//) 4 Hz- 1 according to Refs. l27U37l . 

7 The amplitude of GWB spectrum might be reduced by a factor of 2 ~ 5 according to the recent numerical simulations 1 1 (I . In case 
A, however, the noise contributions in SNR 1311 used for reconstruction analysis are basically determined by the k = components of 

self-correlation signals, i.e., = u CjjfiCjjfi. Hence, a slight change of the GWB amplitude would not alter the final results. 
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FIG. 6: Signal-to-noise ratio for each Fourier component of correlation signals, (S/N)k estimated at / = 0.1 (/ ~ lmHz). The 
histogram depicted in solid line indicates the case with the amplitude of Galactic GWB given by Sl /2 = 5 x 10" 19 Hz- 1/2 at / = 
lmHz (case A), while the histogram in dashed line represents the signal-to-noise ratio for the case with Sj/ 2 = 5 x 10~ 18 Hz~ 1//2 
(case B). The thin-dashed, and thick-dashed lines respectively denote the lines of (S/N)k = 1 and 5. 



AA and EE as well as for the cross-correlation signal AE are generally good compared to the cross-correlation data 
AT and ET. With sufficient higher SNR of (S/N) k > 5, the available Fourier components of AA-, EE- and AE- 
correlations become k = —2, —1, 0, +1, +2 in both realistic and optimistic cases. This is consistent with the previous 
estimates [l8j- 8 On the other hand, with a large amplitude of GWB spectrum (case B), only the k = component 
is accessible in the signal combinations of AT and ET. This is mainly due to the fact that the sensitivity of the 
T-variable is poor at low-frequency and thereby the noise contribution N k becomes {CAAfiS"^} 1 ^ 2 or {Cee&S^} 1 / 2 , 
which is much larger than {Caa.oCtt.o} 1 ^ 2 or {CEEfiCrTfi} 1 ^ 2 ■ 

Thus, the available Fourier components of correlation data used for the reconstruction of the skymap would be 
severely restricted in practice. Under such restricted situation, the deconvolution problem of the linear system 1|12|) 
tends to be under- determined. Nevertheless, as it will be shown below, one can determine the £ = 0, 2 and 4 modes 
of multipole coefficients of the Galactic GWB with sufficiently small errors. In addition, with the k — components 
of AT- and i?T-correlations, the odd moments £ — 1 and 3 can be recovered. 



C. Reconstruction of a skymap 



Keeping the remarks on the SNR estimation in previous subsection in mind, we now proceed to a reconstruction 
analysis and test the validity of perturbative reconstruction method presented in Sec. IIIII For this purpose, in addition 
to the analysis in the under-constrained cases (case A and B) mentioned above, we also consider the over-constrained 
case as an illustrative example. 



The precise numerical values of the SNR for AA- and i?_B-correlations slightly differ fr om Ref. llSl , This is mainly because the Euler 
rotation angles for the LISA orbital motion given in Sec III Al are different from those in Il8l . 
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FIG. 7: Reconstructed skymap from the time-modulation signals in the over- determined case. The results were obtained by 
the method based on harmonic-Fourier representation and are plotted in the Galactic coordinate. The left panel shows the 
leading-order result, where £ = 0, 2 and 4 modes are only reconstructed, while the right panel represents the result including 
all reconstructed multipoles (I < 5). 



1. Over-determined case 



(I) Harmonic- Fourier representation Let us first focus on a very idealistic situation that the noise contributions 
are entirely neglected. In such a case, all the components of self- and cross-correlation data Ck are available to the 
reconstruction analysis. In practice, however, it is sufficient to consider some restricted components among all available 
data. Here, to make a skymap with multipoles of I < 5, we use the k = —2 ~ +2 components of self-correlation 
signals AA and EE, the k — —4 ~ +4 components of cross-correlation signals AT and ET, and k = —8 ~ +8 
components of AE-signal (see Table [HJ|. Even in this case, the linear system (|12|) is still over-determined and the 
least-squares approximation by SVD is potentially powerful to obtain the multipole coefficients of anisotropic GWB. 
The procedure of reconstruction analysis is the same one as presented in Fig0 For the output signals of harmonic- 
Fourier representation, we use the data Ck observed at the frequencies / = 0.05 and 0.15, in addition to the data for 
our interest at / = 0.1. Collecting these multi-frequency data, the perturbative expansion form of the vector c(/) is 
specified up to the third order in / and the coefficients are determined (Appendix Q. 

In Fig0 the reconstructed results of multipole coefficients pi m are converted to the projected intensity distribution 
@ and are shown as Hammer- Aitoff map in Galactic coordinate. Left panel shows the skymap reconstructed from the 
lowest-order signals of self- and cross-correlation data c^ 2 \ which only includes the I = 0, 2 and 4 modes, while right 
panel is the result taking account of the leading-order correction c^ 3 \ Comparing those with the expected skymap 
shown in Fig0J the reconstruction seems almost perfect. In Fig. [SJ numerical values of the reconstructed multipole 
coefficients are compared with the true values listed in Table ITTT1 The agreement between the reconstruction results 
(open circles) and the true values (crosses) is quite good and the fractional errors are well within a few percent except 
for P20 and p§o. A remarkable fact is that the monopole and the quadrupole values can be reproduced reasonably 
well despite the presence of the degeneracy mentioned in Appendix^] This is just an accidental result. The (small) 
discrepancy in the "recovered" P20 can be ascribed to the expression l|B4[l . After all, the least-squares method by 
SVD provides a robust reconstruction method when the linear system (|12|) becomes over-determined. 

(II) Time-series representation The successful reconstruction of the GWB skymap can also be achieved by the 
alternative approach based on the time-series representation JSJ. Following the procedure presented in Sec lIII Bl the 
intensity skymap of the GWB is directly obtained and the results taking account of the lowest-order and the leading- 
order contributions to the antenna pattern functions are shown in left and right panels in Fig|^l respectively. Here, to 
create the discretized data set (12611 . the number of grid and/or mesh was specified as M — 16 in time and N = 17 x 32 
in spherical coordinate. The time-series data of antenna pattern functions were numerically generated based on the 
full analytic expressions given in Sec lll Bl under assuming that the arm-length of the three space crafts are rigidly kept 
fixed. The reconstructed skymap reasonably agrees with Fig[7| as well as Fig0] Although the situation considered 



TABLE II: Components of Ck used for reconstruction analysis based on the harmonic-Fourier representation 
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FIG. 8: Reconstructed values of multipole coefficients \pe m \ m over-determined case. These numerical values are evaluated 
in the ecliptic frame. The open circles represent the reconstruction results, while the crosses mean the true values, which are 
obtained from the spherical harmonic expansion of full-resolution skymap (see Table ITTTI . 




FIG. 9: Reconstructed skymap from the time-modulation signals based on the time-series representation. Here, to perform 
the reconstruction analysis, the number of grid and/or mesh for the discretized data set I2(jll was specified as M = 16 in time 
and N = 17 x 32 in spherical coordinate. The left panel shows the lowest-order result in which only the I = 0, 2 and 4 modes 
are included, while the right panel represents the result taking account of the leading-order correction to the antenna pattern 
functions, which includes the multipoles, I < 5. 



here is very idealistic and thus the results in FigEI should be regarded as just a preliminary one, one expects that the 
methodology based on the time-series representation is potentially powerful even when the rigid adiabatic treatment 
of space craft motion becomes inadequate. To discuss its effectiveness, a further investigation is needed. The details 
of the analysis including the effects of arm-length variation will be presented elsewhere. 



2. Under- determined case 



Turning to focus on the analysis based on the harmonic-Fourier representation, we next consider the under- 
constrained case in which the number of available Fourier components is restricted due to the noises (case A and 
B discussed in Sec lIV B|) . FigllOl shows the reconstructed images of GWB intensity map in Galactic coordinate, free 
from the noises but ] restricting the number of Fourier components according to Tabic [H] The top panel is the 
intensity map obtained from the lowest-order signals c^ 2 \ Since the accessible Fourier components are basically the 
same in the lowest-order analysis, the same results are obtained between case A and B. On the other hand, the 
bottom panels of Fig^l] represent the skymap reconstructed from both c' 2 ' and c^ 3 \ which indicate that the different 
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FIG. 10: Reconstructed skymap from the time-modulation signals using the restricted Fourier components in under-determined 
cases (case A, B). These are all plotted in Galactic coordinate. The upper panel shows the result from lowest-order analysis in 
which only the multipole coefficients of I — 0, 2 and 4 modes are considered. The bottom panels are the intensity map taking 
account of the leading-order correction including the multipoles, £ < 5. Left and right panels respectively show the results 
obtained from the case A and B. Note that the available number of Fourier components was restricted in the reconstruction 
analysis according to Table ITT1 

images of GWB skymap are obtained depending on the available number of Fourier components (or the amplitude of 
GWB spectrum); case A (left) and case B (right). In Fig^] numerical values of the reconstructed multipoles pe m in 
ecliptic coordinate are summarized together with the statistical errors. The statistical errors were roughly estimated 
according to the discussion in Appendix El based on the SNR [Eq. 1311) ]. 

It turns out that the case A fails to reconstruct all the dipole moments, while they can be somehow reproduced 
in the case B. This is because the cross-correlation data AT and ET which include the information about 1=1 
modes were not used in the reconstruction analysis in case A. Although the I = 5 modes of multipole coefficients 
were not reproduced well in both cases, their contributions to the intensity distribution are not large. Hence, the 
reconstructed GWB skymap in case B roughly matches the expected intensity map in Fig0]and the visual impression 
becomes better than case A. This readily implies that large amplitude of the GWB spectrum is required for a correct 
reconstruction of a skymap, in practice. However, it should be emphasized that the present reconstruction technique 
can work well in the under-determined cases. Even in the realistic situation with smaller amplitude of GWB spectrum 
(case A), the reconstructed skymap including the multipoles £ < 5 still shows a disk-like structure, which may provide 
an important clue to discriminate between the Galactic and the extragalactic GWBs. 

V. CONCLUSION AND DISCUSSION 

In this paper, we have presented a perturbative reconstruction method to make a skymap of GWB observed via 
space interferometer. The orbital motion of the detector makes the output signals of GWB time-dependent due to the 
presence of anisotropies of GWB. Since the output signals of GWB are obtained through an all-sky integral of primary 
signals convolving with an antenna pattern function of gravitational-wave detectors, the time dependence of output 
data can be used to reconstruct the luminosity distribution of GWBs under full knowledge of detector's antenna pattern 
functions. Focusing on the low-frequency regime, we have explicitly given a non-parametric reconstruction method 
based on both the harmonic- Fourier and the time-series representation. With a help of low-frequency expansion of 
the antenna pattern functions, the least-squares approximation by SVD enables us to obtain the multipole coefficients 
of GWB or direct intensity map even when the system becomes over-determined or under-determined. For illustrative 
purpose, the reconstruction analysis of the GWB skymap has been demonstrated for the confusion-noise background 
of Galactic binaries around the low-frequency / = lmHz. It then turned out that the space interferometer LISA free 
from the noises is capable of making a skymap of Galactic GWB with angular resolution up to the multipoles, I = 5. 
For more realistic case based on the estimation of signal-to-noise ratios, the system tends to become under-determined 
and the number of Fourier components used in reconstruction analysis would be severely restricted. Nevertheless, the 
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FIG. 11: Reconstructed values of the multipole coefficients \pe m \ in the under-constrained cases. The numerical values are 
evaluated in the ecliptic coordinates. Left (right) panel shows the result in case A (case B). The open circles with error bar 
represent the reconstruction results, while the crosses mean the true values. The error bars represent the statistical error in 
the signal processing. The evaluation of the error is discussed in Appendix El 



resultant skymap still contains the information of the multipoles up to £ < 5, from which one can infer that the main 
sources of GWB come from the Galactic center. 

Although the present paper focuses on the reconstruction of the GWB skymap from the LISA, the methodology 
discussed in Sec lIIIl is quite general and is also applicable to the reconstruction of low- frequency skymap obtained from 
the ground detectors as well as the other space interferometers, provided their antenna pattern functions. Despite a 
wide applicability of the present method, however, accessible multipole moments of low-frequency GWB are generally 
restricted to lower multipoles due to the properties of antenna pattern functions. This would be generally true in 
any gravitational-wave detectors at the frequencies / < /* = c/(2ttL), where L is arm-length of single detector or 
separation between the two detectors. Hence, with the limited angular resolution, only the present methodology may 
not give a powerful constraint on the luminosity distribution of GWBs. In this respect, we need to combine the 
other techniques such as the parametric reconstruction method, in which we specifically assume an explicit functional 
form of the luminosity distribution characterized by the finite number of model parameters and determine them 
through the likelihood analysis. The data analysis strategy to give a tight constraint on the luminosity distribution is 
definitely a very important issue and it must be considered urgently. Nevertheless, we note that detectable multipole 
moments depend practically on the signal-to-noise ratios. Since the signal-to-noise ratios are determined both from 
the amplitude of GWB and detector's intrinsic noises, a further feasibility study is necessary in order to clarify the 
detectable multipole moments correctly. With improved data analysis strategy, it might be even possible that the 
angular resolution of GWB skymap becomes better than that of LISA considered in the present paper. 

Another important issue on the map-making problem is to consider the reconstruction of skymap beyond the low- 
frequency regime, where the antenna pattern functions give a complicated response to the anisotropic GWBs and 
thereby the angular resolution of GWB map can be improved |l9| . Since the low- frequency expansion cannot be used 
there, a new reconstruction technique should be devised to extract the information of anisotropic GWBs. Further, 
in the case of LISA, the effect of arm-length variation becomes important and the rigid adiabatic treatment of the 
detector response cannot be validated. Hence, as emphasized in Sec lIII Bl it would become essential to consider the 
reconstruction method based on the time-series representation, in which no additional assumptions for spacecraft 
configuration and/or motion are needed to compute the antenna pattern functions. The analysis concerning this issue 
will be presented elsewhere. 



18 



Acknowledgments 

We would like to thank N. Seto for valuable comments on the estimation of signal-to-noise ratios. We also thank 
Y. Himemoto and T. Hiramatsu for discussions and comments, T. Takiwaki and K. Yahata for the technical support 
to plot the skymap. The work of H.K. is supported by the Grant-in- Aid for Scientific Research of Japan Society for 
Promotion of Science. 



APPENDIX A: MULTIPOLE COEFFICIENTS OF ANTENNA PATTERN FUNCTIONS IN 

LOW-FREQUENCY REGIME 

In this appendix, based on the analytic expression (21) of paper I, the multipolc moments for antenna pattern 
functions of optimal TDIs are calculated at detector's rest frame. Using the low-frequency approximation / -C 1, we 
present the perturbative expressions up to the forth order in /. 

To evaluate the antenna pattern functions, we must first specify the directional unit vectors a, b and c, which 
connect respective spacecrafts. Here, we specifically choose (Fig. 1 and Eq.(28) of Paper I): 

V3 1 V3 1 

a = -— x+-y, b= y, c = — x + - y, 

where the vectors x and y respectively denote the unit vectors parallel to the x- and j/-axes in detector's rest frame 
(see Fig.l of Paper I). Then, based on the expression Ijl3|l . the antenna pattern functions of the optimal TDIs, Tjj 
(/, J = A,E, T) are analytically computed at detector's rest frame. Their multipole moments become 



U 



aim{f)= d9 d(f> sin 6* Y^ m (0, </)) JF(/, 9, </>), (Al) 



which are expressed as a function of normalized frequency / = ///*■ Here, for definiteness, we write down the explicit 
form of the harmonic functions Yg m : 



rr(M) = Pr{cos e) e **. (A2) 

The analytic expressions for the multipole moments a£ m are generally intractable due to the complicated form of 
the antenna pattern functions. In the low-frequency regime, however, the perturbative treatment regarding / as a 
small expansion parameter is applied to derive an analytic expression of multipole moments. Below, we summarize 
the perturbation results up to the fourth order in /: 



_ 20? ? 2 2110r? 4 _ 4 n !6 u l-iy/3 /15tt H 

T AA : a 00 -—f-^^-f, a20 =-y-/ --^-/ a 22 - ^ 

V " I3y/ir H 13(1 -i\/3) FIT u 1 + FIT ~ 2 13(l + i\/3) [Stt ? 4 

a4 °-105 f + ~924^ / ' a42 ~ 5544 VlO 7 ' ° 44 ~ 6" V7G / + 792 V 14 f 



1 l~Tt~ fi (1— «V3) / 37T ~ 4 1 + i-s/3 / /T 



060 - 1848 V 13 /4 ' a62 ~ 4T52 - V 455 /4 ' ° 64 ~ ~^9^ V 1^2 A (A3) 

for self-correlated antenna pattern, Taa- Note that the multipole moments of Tee are related to those of Taa 
through the relations, af^ = a^ A for m = 0, ±6, ±12, ±18, • • • and af^ = — af^ for m = ±2, ±4, ±8, • • • (see Paper 
I). The multipole moments of cross-correlation signals are 

3 + iV3 [5k H 1 [7n f3 _ 13(\/3 + i) 



Tab: a 22 -^^^j-f , a 33 - - f , a i2 - ——J-f, 



\/3 — i fir ~ 2 13(V3 — i) /57T - 4 1 / 7r ~ 3 \/3 + i / 3n » 4 



° 44 6 \J 70^ 792 V 14 ^ ' fl53 18 V 1155 ^ ' 062 4 7 52 V 455 ^ 

(A4) 

396 V 182 J ' y ' 



_ (l + iVS)^ ~ 3 _ V3 + i /37T H _l + i v / 3 / 7r 

^ at: an - — 168 — f > a22 --^r\lT f ' a3l -^rvi4 / 
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«42 



a 55 



v / 3 + ^ R ~ 4 = 17(73-0 

3168 V 5 ; ' fl44 3 1 68 



1 — iy3 / 37r » 3 



72 



154 



«62 



n/3- 



4752 



f4 

J ) «51 



i4 



l + i-s/3 AT 
55 

3168 V 91 



1008 



; 3 , 
/ 4 . 



(A5) 



The multipolc moments of antenna pattern function Tet are also obtained from those of Tat using the relation, 



-(i/y/S) tan (mn/S) af^, as shown in Paper I. 



..ET 

Finally, we also list the multipole moments of self-correlated antenna pattern Ttt, which are all higher-order 
contribution with C(/ 4 ): 



O40 



1584 



1 



a 60 



11088 V 13 



1 



a 66 — — 777 



48 V 3003 



(A6) 



APPENDIX B: REMARKS ON THE DEGENERACY BETWEEN MONOPOLE AND QUADRUPOLE 
MOMENTS FOR LISA MEASUREMENT OF GWB ANISOTROPY 

The reconstruction method based on the harmonic- Fourier representation l|12(l presented in Sec lIII Al has been first 
discussed by Cornish [TEl llt| . Later, Seto & Cooray 0] considered the reconstruction of a skymap in the low- 
frequency limit using the optimal set of TDI signals A and E. In this case, the accessible multipole moments pi m are 
£ = 0, 2 and 4 (see Table [|J. Seto & Cooray explicitly wrote down the expressions for linear equation (|12|l in the case 
of the self-correlation signals (i.e., AA-, EE-correlations) and showed that the output data with sufficient statistical 
significance are only Co, C±± and C±2- Further, they found that the multipole coefficients poo and P20 cannot be 
separately determined due to the degeneracy associated with a specific combination between the Wigner D matrices 
and the antenna pattern functions. 

Here, we explicitly point out the origin of this degeneracy and discuss its influences on the reconstruction of GWB 
skymap. Using the multipole coefficients of antenna patterns in Appendix^ the lowest-order contribution to the 
linear system (|12|) for AA- and EE-correlations becomes 

= I^F { | *» - Wo - tJ^ P4o - 7^ 4 Ct' P, m } , (Bl) 

6 ™>° = I^sf { \ Pm - 117! P20 - iSo P4 ° + ^ m £ 4 ^ P4m } ■ (B2) 

with 6^ A) = 1 — i\/3 = [6 4 A * 4 ]*. Here, we only show the relevant components of Ck which contains the multipole 
coefficients poo and p2o- The above expressions include the multipole coefficients of £ = 4, which are all determined 

— (2) 

separately from the lowest-order contribution of ^4i?-correlation, C° AE k . Thus, apart from the octupole moments, 



equations (|B1(I and (|B2(1 clearly show the presence of degeneracy between the remaining multipole coefficients, poo 
and P20- 

In a language of the least-squares method by SVD, this degeneracy implies that the sub-system in the matrix 
equation (|24|l containing the coefficients j?oo and j?20 becomes under-determined. In this case, the least-squares 
method by SVD cannot correctly produce the approximate solutions for poo and P20 , although it still provides some 
"approximate" solution. From equations IjBlfl and i|B2f) . the only meaningful equation for poo andp20 is now reduced 
to 

C=-± J =(\p O0 - -t-7= P2o\ , (B3) 



4^5"" 14V5 J 

where the numerical constant C represents a collection of the irrelevant terms of £ = 4 modes, which are separately 
determined from the linear system in Ai?-correlation. If one naively applies the least-squares method to the above 
system, a very tight relation is obtained: 

784007- ~ 280\/5l ~ 

Poo = C, P20 = C. (B4) 

The presence of degenerate coefficients may be a big obstacle in constructing the skymap as well as in determining 
the normalization factor of GWB spectrum. In principle, this degeneracy can be broken when we consider the higher- 
order terms of C(/ 4 ). However, these terms are generally small and irrelevant for the reconstruction analysis in the 
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low-frequency regime. In this sense, the reconstruction of low-frequency skymap is, in nature, incomplete and the 
other additional information for monopolc or quadrupolc moment is required to make a full skymap. Nevertheless, 
it should be emphasized that with a high signal-to-noise ratio, the other remaining multipole coefficients can be all 
determined by the least-squares solution by SVD, independently of the above degeneracy. Moreover, in cases with 
Poo 3* P20, which is usually satisfied, the least-squares solution (|B4|I provides a modest estimate of the degenerate 
coefficients poo and P20 because of the hierarchy of the coefficients in Eq. I|B4(I . This point has been explicitly 
demonstrated in Sec II VI 



APPENDIX C: DATA SETS AND LEAST-SQUARES METHOD 



In this appendix, we describe in more details how to obtain the multipole moments of GWBs from many data 
streams by applying the least-squares method. Here we specifically focus on the situation considered in Sec. IIV C II 
That is, the data sets that we use are (i) the k — —2 ~ +2 components of self-correlation signals AA and EE, (ii) 
k = —8 ~ +8 components of AE-signal, and (iii) the k = —4 ~ +4 components of cross-correlation signals AT and 
ET. According to the general strategy given in Sec. IIII Al (see Eq. JUJ), we first combine the leading data streams 
of (i) and (ii) which correspond to i = 2 in Eq. I|23[l . The combined data sets consist of the following matrices. 



/ c (2) 



+2.AA 
°+l,AA 

rv(2) 



c 



(2) 

+2, EE 



\ c {2) I 



/o * 
0*0 
0*00 

* 
* 
0*0 
0*00 

* 



/ r (2) \ 




c 



+8,AE 
+7,AE 



(2) 

O.AE 



M 2 ) 

^-l.AE 

\ r (2) ' 

\ ^-S.AE 



V 



/ 



( P44 \ 

P43 



P40 



P4.-3 
\ P4.-4 / 



(CI) 



where * represents a non- vanishing complex number. Note that for simplicity we have not included the data Cq aa 
and C^ EE , which contain degeneracy between the multipole moments (see Sec. EJ. The matrices in the right hand 

side correspond to in Eq. I|23() . These sparse forms of matrix are typical for the present problem. We perform 
the singular value decomposition with respect to these sparse matrices, and construct pseudo-inverse matrices which 
give the least-squares solutions of pt m . In a similar way the least-squares solutions of £ = odd modes are obtained 
from the following matrix equation, 



/ r (3) \ 



M 3 ) 
u -8,iB 



C 



(3) 

4, AT 



M3) 
U -4.1T 

Ly 4._ET 



V 



A 
A 



(3) 
AE 

(3) 
AT 

(3) 
ET 



P10 
Pi -1 
P33 



P3.-3 
P55 



J 



(C2) 



V P5-5 I 

where we have symbolically written the matrix A' 3 ' due to space limitation 



I M3) I 



APPENDIX D: MULTIPOLE COEFFICIENTS FOR GRAVITATIONAL- WAVE BACKGROUNDS 

Here, we give the multipole coefficients for normalized intensity distribution P(ft) of galactic GWB presented in 
Sec IIV Al To evaluate this, we first create the intensity map P(fl) with 129 x 256 regular grids of celestial sphere (9, <f). 
Then, the spherical harmonic expansion of the intensity map is numerically carried out using the SPHEREPACK 3.1 
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TABLE III: Multipole coefficients for gaiactic GWB. pi m with m < is given by pe, n = ( — l) m Pi _ 



1} 

t 


m 




lm^ m J 








0.282 





1 





-0.0215 





1 


1 


0.00876 


-0.156 


2 





-0.0681 





2 


1 


-0.0828 


0.0216 


2 


2 


-0.180 


-0.0124 


3 





0.0231 





3 


1 


0.00419 


0.0648 


3 


2 


0.0212 


0.0544 


3 


3 


-0.0170 


0.135 


4 





0.0120 





4 


1 


-0.0205 


-0.0224 


4 


2 


0.0807 


-0.00274 


4 


3 


0.0936 


-0.0198 


4 


4 


0.139 


0.0194 


5 





-0.0172 





5 


1 


0.000230 


-0.0272 


5 


2 


-0.0250 


-0.00640 


5 


3 


-0.00336 


-0.0693 


5 


4 


-0.0177 


-0.0719 


5 


5 


0.0222 


-0.1.13 



package |35| . Table ITTT1 summarizes the multipole coefficients pe m up to I = 5, which are relevant to the analysis in 
Sec|lV| Note also that p e< _ m = {-l) m p* tm . 

To characterize the contribution of the £-th moment to the galactic GWB, we introduce the angular power defined 

by 

E . ( D1 ) 

n=-l ) 

which is rotationally invariant Figure^lshows the normalized angular power, cti/gq up to I — 30. The dominant 
contribution to the intensity of galactic GWB comes from the multipolcs with t < 4, however, the asymptotic behavior 
at higher multipoles is very slow and can be fitted by a e /a - 1.85 e -° 005 7^ (dotted line in FigU21), which turns out 
to be a good approximation even to much higher multipoles, £ ~ 100. 





FIG. 12: Normalized angular power of galactic GWB anisotropy, ae/ao as a function of multipole I. 
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APPENDIX E: ESTIMATION OF STATISTICAL ERROR IN RECONSTRUCTION ANALYSIS 

In the reconstruction analysis based on the harmonic-Fourier representation, the statistical errors plotted in FigfTD 
are roughly estimated as follows. In the presence of the noises, the least-squares solution given in equation l|25|) 
becomes 



rapp 



A«] + -{c«+s«}, (El) 



where the additional term si 1 ^ represents the noise contributions to the i-th order coefficient of perturbative expansion 
for c(/). The root-mean-square amplitude of the error Ap^ is then defined by taking the ensemble average of the 
random noises as ApW = (\p^ — (p^)| 2 ), which gives the errors in multipole coefficient Apg m . The j-th components 
of the vector ApW becomes 



|Ap W .1 

I -t^approx j I 





+ 


A w" 




jk 





, (\&)- (E2) 




Here, the variance (Is^l 2 ) roughly corresponds to the quadrature of the vector cW divided by the signal-to-noise 
ratio: 

r „w V 

(E3) 

In the above expression, the vector (s/r) represents the SNR, each component of which is the quantity (S/N)k defined 
in (|31|) just for the same component of the vector c' 1 '. Notice that the factor a is multiplied in equation i|E3(l by hand. 
The reason why we have introduced the factor a is as follows. First note that the quantity (S/N)k basically reflects 
the signal-to-noise ratio for the most dominant term in the perturbative expansion for the signal Ck(f) (see Ea. (|23f) '). 
For the J 4£'-correlation, (S/N)k gives the signal-to-noise ratio for second-order term, i.e., c^ 2 \ For the AT-correlation, 
(S/N)k reflects the signal-to-noise ratio for c^ 3 ). In our reconstruction analysis, the higher-order contributions to the 
AE-correlation, c^ 3 ) is used to make the skymap with multipoles I < 5, which can be estimated by analyzing the 
multi-frequency data. In general, the signal from higher-order contribution is weaker than the lowest-order term, 
leading to the calibration error. The significance of this error would be enhanced by the factor roughly proportional 
to j^ 4 " 2 -' for the i-th order terms of AE-correlation. Hence, in order to mimic this, the factor a is multiplied and set 
to a = Z^ 1-2 -*. In the case examined in Fig^] we adopt a — f^ 1 = 10 for third-order cross-correlation data C^e k- 
Otherwise we set a = 1. As a result, statistical errors of £ — 3 modes become larger than those in £ =even modes 
(see Fig lllH . Note that the multipole coefficients with £ = 1 in case A and with £ — 5 in both case are completely 
degenerate and cannot be recovered by reconstruction analysis. Hence, the statistical errors were not evaluated. 

APPENDIX F: COMPUTATIONAL METHOD OF MULTIPOLE COEFFICIENTS 

In this Appendix, we give a brief description of the numerics of calculating the multipole coefficients using the 
SPHEREPACK 3.1 package 35]. The traditional spherical harmonic transform of a scalar function is 

= X>miWM). (fi) 



The spherical harmonics are given by equation (|A2I) and the associated Legendre polynomials are 

On the other hand, numerical computation of spherical harmonic expansion is performed with the subroutine shaec 
in SPHEREPACK 3.1 package, which directly gives the following expansion form: 

I — oo oo 

V{9,4>) = J 2 ^2 ^2 ' P ™( cos6 07ftn [a im cos m(f>- /3 £m sinm0] , (F3) 

" 1=0 m=0 
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where the prime notation on the sum indicates that the first term corresponding to m = is multiplied by the factor 
1/2. To relate the coefficients cet m and f3g m with pi m , the expansion I|F1[) is compared with the expression {FJJ), 
leading to 

Pirn = t/-(aim + !'Wi 
{-l) m pi,- m = ^(a em -i(3 (m ), (F4) 

for m > 0, and 

2 a eo, (F5) 

for m = 0. Hence, with a help of these expressions, one can read off the multipole coefficients pi m from the expansion 
formula JF2J). Notice that the associated Legendre polynomials used in the SPHEREPACK 3.1 package differ from 
(|F2|I by a factor (— l) m so that one must multiply this factor to the transformation law l|F4j) and l|F5|l . 
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